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We introduce a numerical method to study random Boolean networks with asynchronous stochas- 
tic update. Each node in the network of states starts with equal occupation probability and this 
probability distribution then evolves to a steady state. Nodes left with finite occupation probability 
determine the attractors and the sizes of their basins. As for synchronous update, the basin entropy 
grows with system size only for critical networks, where the distribution of attractor lengths is a 
power law. We determine analytically the distribution for the number of attractors and basin sizes 
for frozen networks with connectivity K = 1. 
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Networks have proven to be useful in many fields. In 
particular for biological high-throughput experiments to 
construct regulatory networks there are immense chal- 
lenges how to interpret their behavior. Relevant criteria 
for models include conceptual simplicity, computational 
tractability, and robustness to uncertainties in the data. 
Boolean networks are candidates for representing classes 
of behaviors observed in large regulatory networks (see 
e.g. Ref. [l|). Biochemical details like reactions rates and 
concentrations (which are often unknown) are discarded 
for a simpler description in which genes or functional 
sets of genes are "on" or "off"; for reviews see [2|, 
In addition to describing specific small regulatory net- 
works d, Q , Boolean networks have been used to repre- 
sent other complex dynamical systems such as neural Q , 
or evolving networks. 

As canonical examples of disordered systems, ran- 
dom Boolean networks (RBNs) [8] have attracted the 
attention of physicists over four decades 0, RBNs 
are directed graphs consisting of N Boolean elements, 
where each element receives input from K distinct el- 
ements. The value of the j*'' element, Ui, evolves ac- 
cording to a random Boolean function of its K inputs: 
ai{t + 1) = fi{(Ji^{t),ai^{t), . . . ,ai^{t)), where (Ji.{t) is 
the value of the input to element i at time t. Here we 
choose the function fi to be zero or one with equal prob- 
ability. The functions fi are fixed for each realization of 
the RBN. 

In the classical RBN (CRBN), elements evolve simul- 
taneously according to a globally synchronized clock. 
States separate into transient states, which cannot be 
reached more than once under the dynamics, and attrac- 
tor states, which can be reached infinitely often in the 
long time limit. Attractor states form closed loops in 
state space. The state space as a whole divides into non- 



overlapping partitions, or basins of attraction. Transient 
states within a given basin all reach the same attrac- 
tor at long times. By considering small perturbations 
to states within an annealed approximation, Derrida and 
Pomeau found that K = 2 CRBNs are critical: they 
separate ordered {K < 2) and chaotic K > 2 phases, 
where distances between nearby trajectories vanish or 
diverge, respectively. Numerous investigations have fo- 
cused on the lengths and numbers of attractors and their 
basins of attraction 0, In particular, Krawitz 

and Shmulevich [l3| found that the basin entropy, which 
measures the variations in the size of the partitions of 
state space, increases with system size only for critical 
CRBNs. 

However, most real world systems do not evolve ac- 
cording to a globally synchronized clock, leaving seri- 
ous doubts on the applicability of CRBNs 18-20] as 
useful models. Here we consider the opposite extreme 
and study asynchronous, stochastic random Boolean net- 
works (ARBNs), where at each time step, a single ran- 
domly selected element is updated. Previous studies 
have found fundamental differences between CRBN and 
ARBN ensembles. For instance Ref. 2l| argued, us- 
ing analytic methods, that shifting from synchronous to 
asynchronous dynamics drastically reduces the number 
of attractors for connectivity K = 2. 

Here we show, using a novel computational method, 
that certain features of RBNs are the same for both ex- 
tremes and thus are generic with respect to the dynamics. 
In particular the basin entropy increases with system size 
only for critical ARBNs. Hence it is a robust detector of 
critical behavior for both CRBNs and ARBNs. Further, 
the distribution of attractor lengths is a power law for 
K — 2 networks in both cases. Finally, we show that all 
attractors \n a. K = 1 ARBN have the same length, and 
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all their basins have the same size. We use this result 
to derive analytically the distribution of the number of 
attractors and the average basin entropy for ensembles 
of K = 1 ARBNs as a function of system size A^. 

An attractor in a discrete system with non- 
deterministic dynamics, such as an ARBN, is also a sub- 
set of all possible states that can be reached infinitely 
often in the long time limit starting from a random ini- 
tial condition. Mathematically, an attractor is a set of 
states such that starting from any state within it: (a) 
all other states of the attractor can be reached and (b) 
these are the only states that can be reached. States that 
do not belong to an attractor are transient. Since tran- 
sient states may form loops in the state space network 
(SSN), our definition is not the same as in Ref. [2l|. We 
define the normalized size of an attractor's basin to be 
the chance of reaching that attractor starting from a ran- 
domly chosen initial state. Note that unlike CRBNs, an 
ARBN can reach different attractors starting from the 
same initial state, i.e. the transient states in different 
basins can overlap. In addition, attractors do not nec- 
essarily form simple loops. However, different attractors 
are still composed of non-overlapping subsets of states. 

We find attractors and their basins using the following 
method: First we construct the entire SSN oi Af — 2^ 
states for an ARBN of size N, by representing the states 
as nodes. A directed link, with weight points from 

a node to each of its N images reached by choosing, in 
turn, one element for update. An image can be the state 
itself or some other state with Hamming distance one 
away. If some among these N images are identical, we 
add the weights of the corresponding links. Second, our 
algorithm initially assigns the same occupation proba- 
bility p = l/M to each node. Then the algorithm up- 
dates this probability distribution on the SSN in parallel 
by dividing and moving the entire occupation probabil- 
ity on each node along its outgoing links according to 
the weights. The algorithm repeats this step until the 
probability distribution is stationary (23l |. By definition, 
the occupation probability p vanishes for transient states, 
but remains finite for states on attractors. We identify 
different attractors as subsets of these latter states that 
are dynamically connected jl^- The value of p for each 
state is the chance of observing that state in the long 
time limit after starting the ARBN from a random ini- 
tial state. The sum of the densities for all states on an 
attractor is its normalized basin size. 

Before presenting our findings using this method, we 
first discuss the case K = 1, where we derive analyti- 
cal results and compare them with results from the SSN 
method, finding complete agreement. The structure of 
a. K ^ 1 RBN is particularly simple: Boolean elements 
form loops with trees rooted in them. Only four Boolean 
functions exist to update elements: copy, invert, force to 
zero and force to one. Relevant components in CRBNs 
determine the number and lengths of attractors as well 



as their basin sizes @, [111- This is also true for ARBNs. 
Loops of elements with at least one forcing function be- 
come frozen eventually, as do elements on trees rooted 
in frozen loops. Hence, relevant components are loops 
that contain no forcing functions. They are even or odd. 
Even (odd) loops have an even (odd) number of invert 
functions. Their number in ARBN i is rieyenii) {nodd{i))- 

Without loss of generality, consider that all functions 
on even loops are copy functions. Then the two stable 
states are all zeroes (00. ..00) and all ones (II. ..II). Hence 
a single even loop generates two point attractors. By 
symmetry, the two point attractors split the SSN into 
two equal parts. Hence, for an ARBN with rieven loops 
and no odd loops, the number of attractors is 2"'^'" . 
Each has a basin of size 2^^"""=". 

Odd loops are equivalent to a loop with one invert 
function, the rest being copy. Assume that the loop is 
in state So — {0 ... 0} and call the element with the in- 
vert function rj. Sq changes only when 77 is chosen for 
update. Then, the new state Si does not change un- 
til the neighbor on the loop downstream from rj is cho- 
sen. Following this argument, we conclude that a sin- 
gle odd loop of length L in an ARBN has 2L different 
states it can be in. An ARBN with Uodd odd loops 
of lengths Li, L2, ■ ■ ■ ^ Ln^^^, each loop j having Sj un- 
frozen elements downstream of it, has one attractor with 
j-jnodd 21^2"^ states. Thus for ARBN i that contains 
neven{i) and nodd{i) loops, the number of attractors is 
Ai = 2"="""^*). Since the attractor from the odd loops 
enters equally into all basins, each attractor a drains a 
basin of size 6q, — 2^~"='"!"(*). Indeed, we observe numer- 
ically that in a given single ARBN the number of states 
in each attractor is the same for all attractors and that 
all basin sizes are also the same. 

We now use results of Ref. 22] to find the distribution 
of the number of even loops, P{neven), over ensembles of 
ARBNs. These authors derive the probability, (3(m), to 
observe the vector m = {mi,TO2, ■ • ■} relevant loops of 
lengths L = {I, 2, . . .} in a. K = 1 RBN. In our case this 
probability reads 



Q(m) = - 



2 J {N - rhy.N"' 2 



i(I + m/iV)n 



L = l 



where m = X^lLi'^^^- The probability P{n) to 
have n relevant loops (both even and odd) is P(n) = 
Q{m\ X^i'^i — '^)- P{n-even) is obtained from P{n) as 



n \ n 



P(neven) = ^ P{n) 



These expressions do not have closed form solutions, but 
can be evaluated numerically. Fig.[l]shows Pirieven) eval- 
uated from Q{m). Noting that P{neven) = -f'(log2^) = 
P(log2 ^), Fig. m also shows P(log2 A) and P(log2 ^), 
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obtained using the numerical SSN method described 
next, finding agreement. 




FIG. 1: The distribution of the number of even loops 
P{neven), obtained using Q(m) - compared to the distribu- 
tion of the number of attractors, P(log2 A), and the distri- 
bution of inverse basin sizes, P(log2 =^), obtained using our 
SSN method for if = 1 ARBNs with N = 16. The number 
of realizations is 2 x 10'' and the error bars correspond to two 
standard deviations. 



We now consider K >1, and distinguish two measures 
for the length of attractors. The first measure, is the 
actual number of states on attractor a. The second mea- 
sure, Aq, takes into account how often states are visited. 
Defining the stationary occupation probability of state j 
on attractor a as Pj,^, the conditional probability of that 
state given that the dynamics reaches attractor a is 



The attractor entropy, s^, determines Aq, through: 

In Aa = Sq = - ^ Qaj log Qaj ■ 

If all states on attractor a are visited equally often then 
Act = la- If J on the other hand, one state dominates, 
Aa — 5- 1. We observe a strong correlation between and 
Act (data not shown). Since we find clearer - although 
not qualitatively different - behavior using A, we present 
those results below. 

Fig. [2^ shows the probability density function for the 
length of attractors, P(A), for N = 16. It is broad for 
K = 1 and K = 2. For /-C = 6, A is either of order 
unity or of the SSN size TV. This is true for &\\ K > 2 
although subleading effects dominate for 2 < if < 6, and 
no clear separation appears for this TV. A rescaled i-'(A) 
is shown in Fig.[2}3- The cutoff in the length Ac ~ 2^ 
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FIG. 2: (a) P(\) for if = 1, 2, 6 and TV = 16. Both if = 1 and 
if = 2 ARBNs show a broad distribution while all attractor 
lengths for if = 6 are either of order unity or of order 2^. (b) 
Rescaled P(A) vs. rescaled A for if = 2 and iV = 12-22. The 
approximate collapse is consistent with a power law decay up 
to a cutoff which is a stretched exponential in A''. 



the data is consistent with a power law decay P{X) ~ 
A~'^^ with Tx = 1.15 ± 0.05. However, the cutoff changes 
shape for different N over the considered range of system 
sizes so the system is not yet in a regime where subleading 
corrections can be ignored. 

The probability, Pa, to reach attractor a from a ran- 
dom initial state is the sum of all occupational proba- 
bilities. Pa = S'^j, Paj over states on a. Its basin size, 
ba = Mpa ■ Ref. [13 proposed that the average basin en- 
tropy (h) , which gives a summary measure for the varia- 
tion in basin sizes, is a measure of complexity and criti- 
cality in CRBNs. The basin entropy of RBN i is 
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suggested by the dimensional arguments in Ref. 21[ and 



a=l 



Pa In Pa 
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and the average basin entropy over the ensemble is (h) = 
J2f hi where R is the number of reahzations. Fig. [3] 
shows the average basin entropy {h) for ARBNs with 
if = 1, 2, 3, 6 and iV = 8 to 20. This figure shows that 
{h) grows only for K ^ 2 while remaining constant for 
K ^ 2. For K — 1, Fig. [3] also shows {h) derived from 
(5(m), which agrees perfectly with the SSN method. For 



CRBNs [17|, as well, the basin entropy was shown to 
grow with system size only for the critical case K = 2 - 
over the same range of system sizes. 



0.4 
0.35 

A 

0.3 
0.25 
0.2 























— K= 


1 










— K= 


2 










— K= 


4 












— K= 


6 






i — 


— I — 










-i 






— i— 


— 











networks using this method in either case. 

To summarize, we have developed a method to find the 
long-time dynamics of ARBNs. It relies on the flow of 
the occupation probability over the complete network of 
possible states until stationarity is reach. We focussed on 
how the space of states divides into overlapping basins of 
attraction and measured the fluctuations of these fuzzy 
partitions in ensembles of ARBNs. Despite the funda- 
mental difference between ARBNs and classical RBNs, 
we flnd that in both cases the basin entropy increases 
with system size only for critical networks. This suggests 
that "superuniversal" features exist for the fluctuations 
in the structure of state space of RBNs - these are invari- 
ant with respect to the speciflc dynamics. As the dynam- 
ics for real world systems most likely lies in between these 
two extremes, our work suggests that the basin entropy 
may be a relevant and robust signifier of criticality, irre- 
spective of the specific dynamics. Finally, we showed that 
all attractors in an ordered K = 1 ARBN have the same 
length, and that all their basins have the same size. We 
used this to derive analytically the distributions of both 
basin entropy and the number of attractors in if = 1 
ARBNs as a function of system size. 
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FIG. 3: The basin entropy (h) for ARBNs with various K, 
and = 8 - 20. (h) grows with N only for K ^ 2. The 
dashed line shows (h) for K — 1 derived from Q(m). The 
solid black line has a slope of 0.15 and is a guide for the eye. 

Ref. showed that the distribution of attractor 
lengths in CRBNs depends on the sampling scheme. For 
CRBNs with K > 2 the distribution obtained by count- 
ing all attractors is a power law, while that obtained by 
randomly sampling initial states is not. For K — 2 both 
distributions are power laws. This is due to a lack of cor- 
relation between the size of an attractor and the size of 
its basin for K = 2, while for if > 2 the basin size grows 
linearly with attractor length. In the ARBNs studied 
here, the basin size varies little with attractor length for 
all K, so the two distributions are indistinguishable over 
the range of system sizes 8 < A'' < 20. 



As in Ref. [17|, a method based on exact enumeration 
of state space creates a severe restriction on the system 
sizes that can be studied. Indeed for if = 2 CRBNs, 
the mean number of attractors grows faster than any 
power law with system size , but one cannot reach this 
asymptotic result for CRBNs in the range of 8 < < 20. 
The growth we observed for CRBNs is much slower (data 
not shown). For ARBNs, the mean number of attractors 
increases as a power law with system size for if = 2 [2l[ . 
We also do not observe this asymptotic result for the 
range of system sizes studied (data not shown). Hence, 
it is not possible to extract an asymptotic growth law 
for the size dependence of the basin entropy for critical 
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